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Abstract. An incremental evolution equation, i.e. a Master equation in statistical me¬ 
chanics, is introduced for force distributions in polydisperse frictional particle packings. 
As basic ingredients of the Master equation, the conditional probability distributions of 
particle overlaps are determined by molecular dynamics simulations. Interestingly, tails 
of the distributions become much narrower in the case of frictional particles than friction¬ 
less particles, implying that correlations of overlaps are strongly reduced by microscopic 
friction. Comparing different size distributions, we hnd that the tails are wider for the 
wider distribution. 


1 INTRODUCTION 

Quasi-static deformations of soft particles, e.g. glasses, colloids, and granular materi¬ 
als, have been widely investigated because of their importance in industry and science. 
However, their macroscopic behaviors are still not fully understood due to disordered 
conhgurations and complex dynamics [1]. At the microscopic scale, their mechanical 
response is probed as a reconstruction of the force-chain network [2], where non-affine 
displacements of particles cause a complicated restructuring of the network including also 
opening and closing contacts. If a macroscopic quantity, e.g. stress tensor, is dehned as a 
statistical average in the force-chain, its response to quasi-static deformations is governed 
by the change of probability distribution function (PDF) of forces. Therefore, the PDFs 
in soft particle packings have practical importance such that many theoretical studies 
[3l H] have been devoted to determine their functional forms observed in experiments [5] 
and simulations [6]. 

Recently, we have proposed a Master equation for the PDFs as a stochastic approach 
towards a microscopic theory for quasi-static deformations of two-dimensional bidisperse 
frictionless particles [7]. The Master equation can reproduce stochastic evolution of the 
PDFs under isotropic (de)compressions, where the conditional probability distributions 
(CPDs) for the Master equation fully encompass the statistics of microscopic changes 
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of force-chain networks. In addition, any changes of macroscopic qnantities constructed 
from the moments of forces can be predicted by the Master equation. 

In this paper, we generalize our stochastic approach towards wider four-disperse size 
distributions and include friction in the contacts. We determine the CPDs from MD 
simulations and compare them with the cases of polydisperse frictionless and bidisperse 
frictionless particles. First, we explain our MD simulations in Sec. Hi Then, we show our 
results in Sec. Ei and conclude in Sec. 01 

2 METHOD 

We use MD simulations of two-dimensional four-disperse mixtures of frictional soft 
particles. The number of constituents {Ni, N 2 , N^, and W) and the size distribution are 
the same with those used in the experiments of wooden cylinders [8] as listed in Table 
[H where the mass, m, is identical in simulations to be used for the unit of mass. The 
normal force between the particles in contact is given by /“• = k^^Xij — rj^Xij with a normal 
stiffness, /cn, and normal viscosity coefficient, r]^. Here, an overlap between the particles 
i and j is introduced as 

Xij Ri Rj dij ( 1 ) 

with an interparticle distance, djj, and the particles’ radii, Ri and Rj. Thus, the relative 
speed in the normal direction is given by its time derivative, Xij. The tangential force is 
introduced as /f = ktHij — r]tyij which is switched to the sliding friction, = —/r|/“ |, 
when the tangential force exceeds the critical value, i.e. if |/h| > fi\f^j\, where kt = kn/ 2 , 
Vt = w/4, and fi are a tangential stiffness, tangential viscosity coefficient, and friction 
coefficient, respectively. Here, yij and yij are a relative displacement and speed in the 
tangential direction, respectively (see the details in Ref. i)- 

To make static packings of N{= Ni + N 2 + + N 4 ) = 1872 particles, we ran¬ 

domly distribute them in a L x L square periodic box, where no particle touches others 
and the friction coefficient is set to zero, /i = 0 (i.e. we use frictionless particles dur¬ 
ing the preparation of static packings). Then, we rescale every radius as Riit + St) = 
[1 -|- {x — Xra(t)} /X]Ri{t) {i = 1,... ,N), where t, St, x, and x^ait) are time, an incre¬ 
ment of time, a target value of averaged overlap, and the averaged overlap at time t, 
respectively. Here, we use a long length scale A = lO^ct with the mean diameter in 
the hnal state, a, to rescale each radius gently 0. During the rescaling, each radius in¬ 
creases if the averaged overlap is smaller than the target value, Xm(t) < x, and vice 
versa, so that the averaged overlap will hnally converge to x. Note that neither par¬ 
ticle masses nor the ratios between different diameters change during the rescaling, i.e. 
ai(t + St)/aj{t -|- St) = aiit)/aj{t). We stop the rescaling when every acceleration of 
particles drops below a threshold, 10“®/cnd/m, and assume that the system is static. 

^We confirmed that static packings prepared with longer length scales, A = lO^tt and lO^d, give the 
same results concerning their power-law behaviors, i.e. critical scaling, near jamming |10| . while we cannot 
obtain the same results with a shorter length scale, A = lOtt. 
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Figure[T]^a) is a snapshot of our simulation, where the system is static with the averaged 
overlap, x = 3 x 10“®d, and each color corresponds to each constituent as listed in Table 
[U Figure [11(b) shows the complete Delaunay triangulation (DT) of the same packing in 
Fig. UKa), where the red solid lines represent force-chain networks, while the blue solid 
lines connect the nearest neighbors without contacts. In our simulations, the distances 
from jamming are determined by the critical scaling of averaged overlap, x = A{(j) — (pj) 
ra. where the critical amplitude is found to be A ~ 0.25cr as shown in Fig. [2](a). We 
also conhrm other critical scalings of the static pressure divided by the normal stiffness 
(Fig. |2](b)) and the hrst peak value in radial distribution function of scaled distance (Fig. 
|2](c)), where we hnd that the power law dependence on the distance from jamming are 
given by p/k^ = 1.25 x (0 — and gi = 0.31(0 — 0j)“^, respectively @ [TU] . 

Then, we switch the friction coefficient to p = 0.5 and apply an isotropic compression 
to the prepared packings by multiplying every diameter by + 5(f)/(f) such that the area 

fraction increases from 0 to 0 -|- 50. After compression, we relax the system until every 
acceleration of particles drops below the threshold again. 

Table 1: The number of particles per species, W, diameters, (Js^ and masses, mg, of the four kinds of 
particles (s = 1,2, 3,4) resembling values used in experiments with wooden cylinders in Ref. [8]. Each 
color in Fig. [T]are also listed in the last column. 


s 

Ns 

as/a^ 

ms/rrii 

colors 

1 

807 

0.3 

1.0 

green 

2 

434 

0.6 

1.0 

red 

3 

414 

0.8 

1.0 

blue 

4 

217 

1.0 

1.0 

gray 


3 RESULTS 

In this section, we introduce a Master equation for the PDFs of forces as a stochastic 
description of microscopic changes of force-chain networks. First, we study microscopic 
responses of force-chain networks to isotropic compressions (Sec. 13.ip . where we describe 
mean values and fluctuations of overlaps in terms of applied strain increments and dis¬ 
tances from jamming (Sec. 13.2p . We then introduce a Master equation (Sec. 13.3p and 
determine transition rates for the Master equation (Sec. 13.4p . 

3.1 Microscopic response 

At the microscopic scale in soft particle packings, the mechanical response to quasi¬ 
static deformations is probed as reconstruction of force-chain networks, where complicated 

^Here, p/k^ is dimensionless in two-dimension and the scaled distances in radial distribution functions 
are defined as r = / {Ri + Rj). 
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Figure 1: (Color online) (a) A static particle packing with an averaged overlap, x = 3 x where 

the total number of particles is N = 1872. (b) The complete Delaunay triangulation (DT) of the static 
packing in (a): The red solid lines are equivalent to force-chains, where their widths are proportional to 
the magnitudes of contact forces. The blue solid lines connect the nearest neighbors without contacts. 
The gray circles represent the particles. 


particle rearrangements cause the recombination of force-chains including also opening and 
closing contacts. To take into account such changes in structure, we employ the complete 
Delaunay triangulation (DT) of static packings as shown in Fig. [11(b), where not only the 
particles in contacts, but also the nearest neighbors without contacts, i.e. virtual contacts, 
are connected by Delaunay edges. We then generalize the dehnition of “overlaps” from 
Eq. (P to 

X-ij = Ri Rj ^ij 5 (2) 

where Dij is the Delaunay edge length and the overlaps between particles in virtual 
contacts {Ri + Rj < Dtj) are dehned as negative values. Because the DT is unique for 
each packing, contacts and virtual contacts are uniquely determined. 

If we apply an isotropic affine compression to the system, every generalized overlap, 
Eq. (P, (not only contacts, but also virtual contacts) changes to 

xf"" = X,, + 1^6,1, , (3) 

where we neglected the higher order terms proportional to Xij6<p and 50^. However, 
the particles are randomly arranged and the force balance is broken for each particle by 
the affine deformation. Therefore, the system relaxes to a new static state, where non- 
affine displacements of particles cause complex changes of contacts including opening and 
closing contacts. After the relaxation, overlaps change to new values, xh that is 

non-affine responses of overlaps. 
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Figure 2: (Color online) Double logarithmic plots of the (a) averaged overlap, x, (b) static pressure 
scaled by the normal stiffness, p/fenj and (c) first peak of radial distribution function for scaled distance, 
5 i, where the dotted lines are power law fittings, (a) x = O.25(t( 0 — ^j), (b) p/fcn = 1-25 x (</) — 
and (c) gi = 0.31(0 — respectively. 




Figure 3: (Color online) (a) and (b): Sketches of static packings (a) before and (b) after (de)compression, 
where the red and blue solid lines represent contacts and virtual contacts, respectively. The four kinds of 
transitions, (CC) contact-to-contact, (VV) virtual-to-virtual, (CV) contact-to-virtual, and (VC) virtual- 
to-contact, are displayed, (c) A scatter plot of scaled overlaps, where the red and blue dots are and 
^affine p^Q^^ed against respectively, (d) A schematic picture of affine and non-afhne responses of scaled 
overlaps, where the blue and red solid lines represent (in average) and the linear functions, 

respectively. The slope in (CC), Oc, and all the dimensionless lengths, fcg and Vg, are proportional to 7 , 
while the slope in (VV) is negligible, a„ ~ 0. 


As shown in Figs. [3](a) and (b), there are only four kinds of transitions from Xij to xF; A 
positive overlap, Xij > 0, remains as positive, xF > 0, or a negative overlap, Xij < 0, stays 
in negative, x^j < 0, where they do not change their signs and thus contacts are neither 
generated nor broken. We call these changes ^^contact-to-contact (CC)” and “virtual-to- 
virtual (VV)”, respectively. The other cases are that a positive overlap changes to a 
negative one and a negative overlap becomes positive, where existing contacts are broken 
and new contacts are generated, respectively. We call these changes opening and closing 
contacts, or in analogy to the above, “contact-to-virtual (CV)” and “virtual-to-contact 
(VC)”, respectively. 

In the following, we scale the generalized overlaps by the averaged overlap before com¬ 
pression such that scaled overlaps before compression, after affine deformation, and after 
relaxation are introduced as .^ = Xij/x, = xf^'^^/x, and = xF/x, respectively 
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(we omit the subscript, ij, after the scaling). From Eq. ([2D and the critical scaling, 
X = A{(j) — 0j), the scaled overlap after affine deformation is found to be 

= ( 4 ) 

which is a linear function of f with an offset, i?a 7 = {Dij/2A(j))'y, proportional to the 
scaled strain increment, 7 = 5(j)/ (0 — (pj). However, the scaled overlap after non-affine de¬ 
formation, f \ fluctuates around mean due to complicated particle rearrangements during 
the relaxation. 

3.2 Mean and fluctuation 

To describe scaled overlaps after non-affine deformation, we measure their mean 
and fluctuations through scatter plots. Figure (HDc) displays the scatter plot, where the 
four different transitions under compression are mapped onto four regions: (CC) f > 0, 
(VV) < 0, (CV) .^ > 0, < 0, and (VC) .^ < 0, > 0, respectively. In this figure, 

overlaps after affine deformation are described by the deterministic equation (|3D, while 
those after non-affine deformation distribute around mean values with finite fluctuations. 
The differences between affine and non-affine responses are always present, but not visible 
if the applied strain is small or the system is far from jamming, i.e. if 7 1 , while f 

deviates more from and data points are more dispersed if 7 ^ 1 . 

In the same way as affine response, Eq. (|3D, we describe the mean values of in (CC) 
and (VV) regions by linear functions of 

msif) = {as + l)^ + bs , (5) 

where the subscripts, s = c and v, represent the mean values in (CC) and (VV), re¬ 
spectively. We also introduce standard deviations of from their means as Vg (which 
are almost independent of ^). Then, the systematic deviation from the affine response is 
quantified by the coefficients, a*, bg, and Vg, as summarized in Fig. [3](d). Note that the 
affine response, Eq. (|3D, is obtained if = Us = 0 and bg = Ba'y. Except for ~ 0, all 
the coefficients linearly increase with the scaled strain increment (Fig. HD, where all data 
with a wide variety of df and f — (j)j collapse onto linear scalings, 

Og = Ag'y , bg = Bg-i , Vg = Vg'j , ( 6 ) 

between the fitting range, 10“® < 7 < 5 x 10“^, with the scaling amplitudes, Ag, Bg, 
and V, listed in Tabled! We observe that ~ 0 and ~ Ba{— 1.3 in average) such 
that virtual contacts almost behave affine in average except for their huge fluctuations 
(W S> Vc). In contrast. Be is always smaller than Ba such that mc{f) intersects 
at = {Ba — Be) /Ac — 1.3, which is independent of 7 . This leads to small responses, 
< ^affine, of Small Overlaps, ^ and vice versa, implying preferred tangential and 

hindered normal displacements as a sign of non-affine deformations m 
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Figure 4: (Color online) Double logarithmic plots of the coefficients for mean values and fluctuations of 
scaled overlaps as functions of the scaled strain increment, (a) Oc, (b) (c) Vc, (e) by, and (f) Vy, and a 

semi-logarithmic plot of (d) a„, where we apply strain increments 6(j) = Ax 10“^, 4 x 10“®, and 4 x 10“® 
(as indicated by the different symbols in the legend of (a)) to the static packings with (l) — (j)j = 1.2 x 10“^, 
4 X 10“^, 1.2 X 10“^, 4 X 10“^, 1.2 x 10“^, and 4 x 10“"^. The solid lines represent the linear scalings, 
Eq. ®, between the fitting range, 10 ® < 7 < 5 x 10 


In contrast to (CC) and (VV), the data of in (VC) and (CV) are concentrated in 
narrow regions (the inside of the dashed lines in Fig. |3](d)), whereas linearly increases 
with ^ in (VC) and there is no data of in (CV), i.e. affine responses do not generate 
opening contacts. 


Table 2; Scaling amplitudes in Eq. ([6]), Ag, Bg, and Vg, for polydisperse frictional particles. The g-indices 
for the CPDs in (CC) and (VV), where and q™ are the q-indices for bidisperse frictionless 

particles [7], polydisperse frictionless particles, and polydisperse frictional particles, respectively. 


c A R T/ r.BL „PL „PN 

c 0.39 0.81 0.30 1.13 1.72 1.19 

V 0.00 1.32 1.22 1.39 1.96 1.09 


3.3 Master equation 

The reconstruction of force-chain networks attributed to the changes, (CC), (VV), 
(CV), and (VC), is well captured by the PDFs of scaled overlaps. Here, we introduce 
the PDF as with the subscript, 0, representing the area fraction in the system. 

Because the total number of Delaunay edges is conserved during deformations, the PDFs 
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Figure 5: (Color online) The PDFs of scaled overlaps before compression, P^{^) (the squares), after 
affine deformation, (the triangles), and after non-affine deformation, (the circles), 

where the inset is a magnification of the PDFs of negative scaled overlaps. 


are normalized as = 1 - 

In our previous study of bidisperse frictionless particles [7], we found that the affine 
deformation, Eq. (0]), just shifted the PDF to the positive direction, while the non-affine 
deformation broadened the PDF in positive overlaps and generated a discontinuous “gap” 
around zero [§. We also find that the PDF of negative overlaps after non-affine deforma¬ 
tion is comparable with that after affine deformation. In our simulations of polydisperse 
frictional particles (/i = 0.5), we observe similar results to our previous study. Figure [5] 
displays the PDFs for polydisperse frictional particles, where the PDF before compres¬ 
sion, has a discontinuous gap around zero. The affine deformation pushes the 

PDF towards the positive direction as , where the discontinuity is smoothed 

out because of the polydispersity in the system. After non-affine deformation, the PDF 
widens in positive overlaps as P(^+ 50 (^'), while there is no significant difference between 
P<j>+S 4 >{^‘^^^'^) and P^+S(i>{0 in negative overlaps (the inset in Fig. [5]). 

To describe such the non-affine evolution of the PDFs, we connect the PDF after non- 
affine deformation to that before compression through the Chapman-Kolmogorov equation 

ra. 

/ oo 

Wi^'l^P^d^ , (7) 

■OO 

assuming that transitions between overlaps (from ^ to ^') can be regarded as Markov pro¬ 
cesses. On the right-hand-side of the Chapman-Kolmogorov equation ([7]), the conditional 
probability distribution (CPD) of scaled overlaps, which were ^ before compression, 
is introduced as 1F(^'|^). By definition, the CPD is normalized as lF(^'|^)d^' = 1. 
Then, the Master equation for the PDFs is derived from Eq. ([7]) as [13] 

^Such a discontinuity is specific to static packings, where a corresponding gap has been observed in a 
radial distribution function of glass with zero-temperature |12| . 
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with the transition rate dehned as = lim^^^o hh(.^'|.^)/(50. The hrst and second 

terms on the right-hand-side of the Master equation (0) represent the gain and loss of new 
overlaps, respectively. Therefore, the transition rates or the CPDs fully determine the 
statistics of microscopic responses of force-chain networks. 


3.4 Conditional probability distributions 

We determine the CPDs of scaled overlaps as the distributions of around their mean 
values, rrisi^)- For example, the CPD for affine deformation is given by a delta function, 
fFaffine(f 10 = ^(0-0®"®), which just shifts the PDF by 5 „ 7 , i.e. P4>+s4>iO = 

However, non-affine deformations generate fluctuations of scaled overlaps around their 
mean values so that the CPDs must have hnite widths. In the following, we determine the 
CPDs for polydisperse frictionless particles and those for polydisperse frictional particles 
to examine the effect of particle friction (characterized by the friction coefficient, /i) on the 
statistics of microscopic responses of force-chain networks. We also compare our results 
with our previous work of bidisperse frictionless particles [7] to study the influence of size 
distributions. 

Figures | 6 ](a) and (c) show the CPDs for polydisperse frictionless particles (/i = 0) in 
(CC) and (VV), respectively. As can be seen, all data are symmetric around the mean 
values, and are well collapsed if we multiply the CPDs and the distances from 

the mean values dehned as Eg = f — ms{f) by 7 and I/ 7 , respectively. In these hgures, 
the solid lines are given by 7W^cc(f 10 = fc{^c/l) and 71^^(010 = fvi^v/l) with the 
q-Gaussian distribution [14] . 


fs{x) 


1 

c{qs) 



n{qs)Vg\ 


1 

1—gs 


(9) 


where the functions are dehned as n(t) = (f —3)/(l —t) and c(t) = Vs\/n(t)B (l/2,n(t)/2) 
with the beta function, B(x, y). In Table [2l we list the q-indices which characterize shapes 
of the CPDs. Note that the g-index must be in the range between 1 < g., < 3, where 
the normal (Gaussian) distribution corresponds to the limit, g ^ 1. Figures [ 6 ](b) and 
(d) are semi-logarithmic plots of Fig. |Hl(a) and (c), respectively, where the dotted lines 
are the CPDs obtained from our previous study on bidisperse frictionless particles [7]. 
From these results, we observe that the basic properties of CPDs, i.e. their symmetry and 
self-similarity for diherent 7 , are not ahected by the size distribution. However, we hnd 
that the CPDs for polydisperse particles have much wider tails than those for bidisperse 
particles, i.e. qf^ gf^ in Tabled implying that the spatial correlation of scaled overlaps 
increases with the increase of polydispersity [T5] . 

Figures m a) and (c) display the CPDs for polydisperse frictional particles (/i = 0.5 
with the same polydispersity as in Fig. | 6 ]) for (CC) and (VV) cases, respectively. Note 
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Figure 6 : (Color online) The CPDs for polydisperse frictionless particles multiplied by 7 in (a) (CC) 
and (c) (VV) plotted against scaled distances from mean values, (s = c,v). (b) and (d) are semi- 
logarithmic plots of (a) and (c), respectively. The solid lines are the q-Gaussian distributions with the 
g-indices, (b) = 1.19 and (d) = 1.09, respectively. The dotted lines are the CPDs obtained from 

our previous study of two-dimensional bidisperse frictionless particles [7], where the g-indices are given 
by (b) g®^ = 1.13 and (d) g®^ = 1.39, respectively. 

that all data are symmetric around their mean values and well collapse after the same 
scaling as in Fig. [ 6 l Figures [71)b) and (d) are semi-logarithmic plots of Fig. [71)a) and 
(c), respectively, where the dotted lines are the CPDs for bidisperse frictionless particles 
[7]. Comparing the results in Figs. [H] and [TJ we confirm that the basic properties, i.e. 
symmetry and self-similarity, of the CPDs do not change by particle friction, while the 
tails of CPDs for frictional particles are much narrower than those for frictionless particles, 
i.e. q™ <C gf'" in Table [2l Therefore, the microscopic friction at the contact drastically 
decrease the spatial correlations of scaled overlaps. 

4 SUMMARY 

We have introduced a Master equation for the PDFs of forces in two-dimensional 
polydisperse frictionless and frictional particle packings. We find that the microscopic 
response to isotropic compression increasingly deviates from the affine prediction with 
the increase of the scaled strain increment, 7 , in agreement with our previous study of 
bidisperse frictionless particles [7]. The scaling amplitudes for mean values of scaled 
overlaps are almost the same as those for bidisperse frictionless particles [7], implying 
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Figure 7: (Color online) The CPDs for polydisperse frictional particles multiplied by 7 in (a) (CC) 
and (c) (VV) plotted against scaled distances from mean values, (s = c, u). (b) and (d) are semi- 
logarithmic plots of (a) and (c), respectively. The solid lines are the q-Gaussian distributions with the 
g-indices, g™ = 1.19 and g™ = 1.09, respectively. The dotted lines are the CPDs obtained from our 
previous study of two-dimensional bidisperse frictionless particles [7], where the g-indices are given by 
(b) g®^ = 1.13 and (d) g®^ = 1.39, respectively. 

that the degree of non-affinity depends on neither particle friction nor size distribntion 
in average. In addition, symmetry and self-similarity of the CPDs, which govern the 
statistics of microscopic changes of the force-chain network, are not affected by friction 
and size distribntion. In contrast, the shapes and widths of the CPDs are different for 
different systems: By nsing the g-indices to describe their shapes in all cases, we observe 
that the g-index for polydisperse particles is larger than bidisperse particles, which means 
that the spatial correlation of scaled overlaps increases with polydispersity (the g-index 
not eqnal to one means a non-Ganssian distribntion). On the other hand, the g-index 
becomes nearly eqnal to one if microscopic friction is considered, implying a dramatic 
decrease of spatial correlation. 

In conclusion, the basic properties of the Master equation, i.e. linear scalings of mean 
values and fluctuations of scaled overlaps, and the symmetry and self-similarity of the 
CPDs, are insensitive to particle friction and size distributions, in the cases tested in this 
paper, while the shape of CPDs depends on both friction coefficient and polydispersity, 
implying different spatial correlations of scaled overlaps for different material properties. 
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